library(tidyverse)
library(here)
library(janitor)
library(raster)
library(sf)
library(tmap)
library(tmaptools)
library(gstat) # variogram
Grand Canyon GeoTIFF from USGS: https://pubs.usgs.gov/ds/121/grand/grand.html
# Read in the data with raster::raster()
gc_dem <- raster(here("data","gc_dem.tif"))
# Look at it with base plot()
plot(gc_dem)
# Check CRS & bounds:
gc_dem@crs # Shows CRS: WGS84
## CRS arguments:
## +proj=utm +zone=12 +datum=WGS84 +units=m +no_defs +ellps=WGS84
## +towgs84=0,0,0
gc_dem@extent # Shows extent (bounds)...notice that these seem odd (not units)
## class : Extent
## xmin : 229520
## xmax : 412580
## ymin : 3982160
## ymax : 4100630
wgs84 = "+proj=longlat +datum=WGS84 +ellps=WGS84 +no_defs" # Just have this ready to copy/paste
# Reproject
gc_reproj = projectRaster(gc_dem, crs = wgs84, method = "bilinear")
# Then check: aha, now degrees we're used to
gc_reproj@extent
## class : Extent
## xmin : -114.0417
## xmax : -111.9681
## ymin : 35.94488
## ymax : 37.04945
bounds <- as(extent(-112.4, -112.0, 36.1, 36.3), 'SpatialPolygons')
# Make the projection for "bounds" the same as for "gc_reproj":
crs(bounds) <- crs(gc_reproj)
# Then crop gc_reproj by the new bounds polygon:
gc_crop <- crop(gc_reproj, bounds)
# Look at it:
plot(gc_crop)
Want to resample? Use raster::aggregate() to create lower res (larger cell) rasters.
See ?aggregrate (default is mean, fact is number of cells in each direction; can set 2 if wanting x/y to differ for aggregation)
# Aggregate:
gc_agg <- aggregate(gc_crop, fact = 10)
# Then look at it:
plot(gc_agg)
# First, convert to a data frame:
gc_df <- as.data.frame(gc_agg, xy = TRUE)
# That `xy = TRUE` is important: retains the lat/lon information!
# View(gc_df)
ggplot(data = gc_df, aes(x = x, y = y)) +
geom_raster(aes(fill = gc_dem)) +
coord_quickmap() +
theme_minimal() +
scale_fill_gradientn(colors = c("purple",
"magenta",
"orange",
"yellow",
"white")
)
Let’s say we know that in this region, a certain species will only grow between 1000 and 1500 ft elevation. Create a subset of gc_crop that includes that habitat:
# First, make a copy
gc_hab <- gc_crop
# Set any cells outside of 1000 - 1500 to NA
gc_hab[gc_hab > 1500 | gc_hab < 1000] <- NA
# Plot.
plot(gc_hab)
# Cool!
Let’s make a bit nicer map with tmap:
tmap_mode("view") # Set to interactive viewing
# Make tmap:
tm_shape(gc_hab) +
tm_raster(legend.show = FALSE, palette = "plasma")
# In console, run tmaptools::palette_explorer() to view a Shiny app with other palettes!
read_sf:ks_counties <- read_sf(here("data",
"ks_counties",
"ks_counties_shapefile.shp"))
# View(ks_counties)
# Base plot
plot(ks_counties)
# Check CRS:
st_crs(ks_counties) # hmmmm none...guess we should set one!
## Coordinate Reference System: NA
# Set to EPSG 4326 (WGS84 datum):
st_crs(ks_counties) <- 4326
# Now check again:
st_crs(ks_counties)
## Coordinate Reference System:
## EPSG: 4326
## proj4string: "+proj=longlat +datum=WGS84 +no_defs"
# And replot:
plot(ks_counties)
# That looks more like Kansas.
ks_rain <- read_csv(here("data","ks_rain.csv")) %>%
clean_names()
But currently, R has no idea that these are spatial points. We’ll convert it using sf::st_as_sf():
ks_sf <- st_as_sf(ks_rain, coords = c("lon", "lat"),
crs = 4326)
# View(ks_sf) (rainfall = amt)
plot(ks_sf)
# Or in ggplot:
ggplot() +
geom_sf(data = ks_counties) +
geom_sf(data = ks_sf, aes(color = amt,
size = amt),
show.legend = FALSE) +
coord_sf() +
scale_color_gradient(low = "yellow", high = "red") +
theme_void()
Now let’s try to predict rainfall all over Kansas based on those rainfall values.
For kriging, we’re gonna switch over to sp functions. Which means we need to get our data to ‘Spatial’ format, instead of sf. Eventually, I think these will all be merged…
ks_sp <- as_Spatial(ks_sf)
# class(ks_sp)
# bbox(ks_sp) to check bounding box of the spatial points
lat <- seq(37, 40, length.out = 200)
long <- seq(-94.6,-102, length.out = 200)
# Then make it into a grid:
grid <- expand.grid(lon = long, lat = lat)
grid_sf <- st_as_sf(grid, coords = c("lon","lat"), crs = 4326)
grid_sp <- as_Spatial(grid_sf)
# Check out your amazing kriging-ready grid:
plot(grid_sp)
# Create the variogram:
ks_vgm <- variogram(amt ~ 1, data = ks_sp)
# Look at it:
plot(ks_vgm)
# Fit the variogram model using reasonable estimates for nugget, sill and range:
ks_vgm_fit <- fit.variogram(ks_vgm, model = vgm(nugget = 0.2, psill = 0.8, model = "Sph", range = 200))
# Plot them both together
plot(ks_vgm, ks_vgm_fit) # Cool! So what are the values
# Just FYI: there are other models (Gaussian, Exponential) - how do those line up?
ks_vgm_gau <- fit.variogram(ks_vgm, model = vgm(nugget = 0.2, psill = 0.8, model = "Gau", range = 200))
plot(ks_vgm, ks_vgm_gau)
# You can check the sum of squares of residuals for each:
attr(ks_vgm_fit, 'SSErr') # 0.00214 (and could compare to other models...)
## [1] 0.002138374
# We'll stick with the Spherical model:
ks_vgm_fit # Nugget = 0.102, sill = 0.954, range = 235
## model psill range
## 1 Nug 0.1021677 0.0000
## 2 Sph 0.9537363 235.1416
Now, kriging!
ks_krige <- krige(amt ~ 1, ks_sp, grid_sp, model=ks_vgm_fit)
## [using ordinary kriging]
spplot(ks_krige, "var1.pred")
Let’s get it back into a format we’re used to (data frames & sf objects):
# Make a data frame from kriged predictions:
ks_df <- data.frame(ks_krige@data["var1.pred"], ks_krige@data["var1.var"],
ks_krige@coords) %>%
rename(longitude = coords.x1,
latitude = coords.x2)
# Convert to sf object:
rain_sf <- st_as_sf(ks_df, coords = c("longitude","latitude"))
st_crs(rain_sf) <- 4326
# Get Kansas outline to crop:
ks <- read_sf(dsn = here("data","states"),
layer = "cb_2017_us_state_20m") %>%
dplyr::select(NAME) %>%
filter(NAME == "Kansas") %>%
st_transform(crs = 4326)
plot(ks)
# Find the intersection of the two:
rain_sf_ks <- st_intersection(rain_sf, ks)
ggplot(rain_sf_ks) +
geom_sf(aes(color = var1.pred)) +
scale_color_gradientn(colors = c("white","yellow","magenta","purple")) +
theme_minimal()